  function [mhm_mat,rej_mat,weightmat]=modhmean_new(vecpar, vecprod, alph); 
%  
% function [mhm,rejrate,weightmat]=modhmean_new(vecpar, vecprod, alph); 
%  Modified Harmonic Mean 
%  Computes the modified harmonic mean as described in 
%  Fernandez-Villaverde and Rubio-Ramirez
%
% Inputs
% -----
%_ vecpar matrix of simulation parameters 
%   [nsim x number of parameters]
%   each row corresponds to a draw of the parameters with MCMC methods 
%_ vecprod matrix containing sum of LOGlikelihood and LOGprior
%   [nsim x 1]
%_ alph px1 vector of chisquared rejection rates 
%   
% Output
% ------
% mhm       p x 2 vector of Log Modified Harmonic Means for the values of alph in the second column
% REJRATE   P x 1 vector with the rejection rates for each value of alph 
% Replaces modhm.m and modhmean.m 
%
% Alejandro Justiniano 10/31/05 
% =================================================================

if nargin < 3; 
    alph = [0.05:0.05:0.95]; 
end 
alph=alph(:) ; 

[nsim,nump]=size(vecpar); 
if size(vecprod,1)~=nsim; 
    error('Dimensions do not match'); 
end; 

% ============================
% Moments of the draws Normal 
% ============================
meanp  =mean(vecpar); 
meanmat=repmat(meanp,[nsim 1]); 

% Deviations from the mean [nsim x npar] 
devm=vecpar-meanmat; 
clear meanmat; 

% Sample variance covariance matrix 
vars=(devm')*devm/nsim; 
varin=vars\eye(nump);
%varin = inv(vars); 
vardet = mydet( vars );  

% Use the determinant of the original covariance matrix 
% Hence the minus sign 
cons_pdf = -0.5*(nump*log(2*pi) + vardet); 

% ===========================================
% Vector of deviations and Normal pdf (log)
% ==========================================
chpr=zeros(nsim,1);
logpdfN = zeros(nsim,1); 
ii=1; 
for ii=1:nsim; 
    chpr(ii)=devm(ii,:)*varin*(devm(ii,:)');  
    logpdfN(ii) = cons_pdf - 0.5*chpr(ii); 
end; 

% Storage Matrices 
% =================
nalph=length(alph); 
mhm_mat = [zeros( nalph, 1) alph]; 
rej_mat = zeros(nalph,1);  
[junk,minalph]= min( alph ); 
clear junk; 

cval=chi2inv(alph,nump); 

difrat = logpdfN - vecprod ; 
temp = sort( difrat ); 
temp_alph = round( length(difrat)*max(alph) ); 
cc    = mean( temp(1:temp_alph) ); 
clear temp temp_alph; 
%cc = max( vecprod );
cc_mat = cc*ones( nsim , 1); 

ii=1; 
for ii = 1:nalph;     
    
    index=find( chpr <= cval(ii) );
    nuse=length(index);     
    nac=nuse/nsim; 
    rej_mat(ii) = 1-nac; 
    marg =  difrat(index) - cc_mat(index) ; 
    marg=exp(marg)/alph(ii);
    if ii == minalph 
        weightmat = marg ; 
    end 
    marg=sum(marg)/nsim  ; 
    
    mhm_mat(ii)    =  -1*(cc + log(marg)); 
   clear nuse nac index marg;  
    
end; 